Multiresolution Wavelet Decomposition I m e Merger of Landsat Thematic Mapper and SPOT Panchromatic Data

Spatially registered Landsat Thematic Mapper (TM) and SPOT (Systeme Pour I'Observation de la Terre) panchromatic images were merged by combining multiresolution wavelet decomposition components from each, and then reconstructing a merged image using the inverse wavelet transform. Three wavelet merging techniques were compared to the intensityhue-saturation merging technique. The comparison results show the wavelet merger providing greater flexibility and the potential for higher accuracy for combining and preserving spectral-spatial information for remotely sensed data and their applications. Introduction The number of commercially available satellite and airborne sensors and the data they are continually increasing. Each sensor has its mission and applications. For many applications, the combined information from multiple sensors provides more comprehensive information by collecting a wide diversity of sensed wavelengths, spatial resolutions, and look-angles. Multiple-sensor exploitation also means manipulation of multiple data sets, some being quite large. This potentially massive data set demands exploration of some means to integrate the relevant information into a more concise, manageable data set. In this paper, we address "fusing" information by exploring a new technique of combining sensor information by using image merging. Image merging refers to image processing techniques that combine two image sets from two or more sensors, forming an enhanced final image. Past merger research mainly explored combining spectral information from a low spatial resolution radiometer with the high spatial resolution information from a wide-band optical sensor. The data combination was performed using simple overlaying, the intensityhue-saturation (IHS) transform merger (Hayden et al., 1982; Carper et al., 1990), component substitution (Shettigara, 1992), and numerous other approaches (Schowengerdt, 1980; Cliche et al., 1985; Tom et al., 1985; Chavez, 1986; Pradines, 1986; Price, 1987; Moran, 1990). These image merging approaches combine the spatial/spectral information from two sensors into one data set, an image. This type of image can be useful in enhancing image mensuration as well as localizing phenomena. If the image merging technique preserves the spectral information, it can be used in spectral classification. Recently, Yocky (1995) tested a novel approach to image merging using multiresolution wavelet decomposition which employs the discrete two-dimensional wavelet transform. He presented the structure and proposed numerous possibilities Systems Analysis Department 111, Sandia National Laboratories, P.O. Box 5800, Albuquerque, NM 87185-0573. in using wavelet decomposition pyramids in image merging. From that research, the "standard" wavelet merger technique was found to outperform the I H ~ merger in preserving the original spectral content while providing comparable spatial resolution. In this paper, we apply discrete two-dimensional wavelet transform image merging techniques to combine Landsat TM data and SPOT panchromatic data. First, we present a brief review of multiresolution wavelet decomposition and wavelet image merging. The "standard" TMISPOT wavelet merge is then presented and compared to the IHS merging technique. We also introduce new algorithms called "additive" and "selective resolution" wavelet mergers. These new wavelet techniques are also compared to the I H ~ merging algorithm. Review of Multiresolution Wavelet Decomposition Multiresolution decomposition (Levine, 1985; Burt, 1989; Wechsler, 1990) provides a simple hierarchical framework for integrating image information. In a pyramidal fashion, image manipulation and analysis can be performed at coarse resolutions proceeding to fine resolutions or vice versa. Mallat (1989) showed how the wavelet transform provides this type of decomposition. The multiresolution wavelet decomposition (MWD) transform is an intermediate representation between Fourier and spatial representations, and it can provide good localization properties in both the spatial and Fourier domains (Daubechies, 1988; Mallat, 1989; Daubechies, 1990). The MWD is computed with a pyramidal algorithm and decomposes a given signal or image into a set of frequency channels of constant bandwidth on a logarithmic scale. The MwD process is presented below starting with the one-dimensional wavelet transform. Daubechies (1988), Mallat (1989), and Chui (1992) provide further mathematical details on wavelets. OneDimensional Wavelet Transform In the following, Z and R denote the set of integers and real numbers, respectively. LZ(R) denotes the vector space of measurable, square-integrable one-dimensional functions f(x) . The multiresolution wavelet decomposition is an increasing sequence of closed subspaces {V,) j E Z which approximate LZ(R). In decomposing the signal f (x) , the resolutions are reduced by a factor of 2 for each level by using a scaling function +(x). The difference in signal at resolutions 21" and 21 Photogrammetric Engineering & Remote Sensing, Vol. 62, No. 9, September 1996, pp. 1067-1074. 0099-1112/96/6209-1067$3.00/0 O 1996 American Society for Photogrammetry and Remote Sensing PE&RS September 1996 (c) (4 Plate 1. Image data set for image merging using part of Sandia National Laboratories, Kirtland Air Force Base, and Southeast Albuquerque, New Mexico. (a) Landsat TM. (b) SPOT panchromatic (copyright, O 1993 CNES). (c) IHS merged image. (d) MWD merged image. can be extracted on a wavelet orthonormal basis of LZ(R). The resolution change is obtained by the first inner product wavelet function is given by Nx), and the wavelet represenwhich acts as a low pass filter: i.e., tation is the orthogonal complement of the original signal space, V,,, denoted as 0 , ) . The decomposition is a new signal h (n) = (4, 4 (x n)) (2) approximation and a detail signal. The signal approximation is given by and by subsampling by two. Using Equation 2 , Equation 1 becomes fz (x) = I: h (2x k) f2,+l (k) where k E Z and (a, b) is the inner product of a and b. The where h(n) = h(-n) . September 1996 PE&RS (c) (dl Plate 2. Image data set for image merging using Albuquerque's south Rio Grande valley region. (a) Landsat TM. (b) SPOT panchromatic (copyright, O 1993 CNES). (c) IHS merged image. (d) MWD merged image. Similarly, the detail signal from the orthogonal projection of f(x) onto O*I is given by d , ( X I = F K q 2 1 ( X I , d (X (k 2n)))(f(x), d,,+, (x 2-1-I k))) . (4) The detail difference between resolutions is obtained by the first inner product which acts as a high pass filter g (n) = (q2-l ( X I , 4 ( x n)) ( 5 ) where g(n) is the quadrature mirror filter of h(n). Using Equation 5 , Equation 4 becomes PE&RS September 1996 where g(n) = g(-n). Thus, the signal approximation at the next lower resolution ( j + 1 + j) is decomposed into a low pass approximation and a high pass detail signal (wavelet coefficients). Perfect reconstruction of the original signal f,+, (x) , requires that the filters h(n) and g(n) have regularity constraints (Daubechies, 1988). The reconstruction is the inverse wavelet transform that takes the form of A,+. (4 = (2k XI fa (kl + g (Zk x) dfd (k) (7) TABLE 1. MWD AND IHS MERGER COMPARISON USING THE AVERAGE ABSOLUTE Although two spatially different images can be transGRAY LEVEL ERROR I N N I N E SEPARATE AREAS. formed to two separate resolution levels in their respective Absolute Gray Level Error pyramids (Step 3, i.e., ' 1 2 for one and '1" for the other), to merge them (Step 4), the image approximation must be the Image Red Green Blue Average correct size to insert into the other sensor's decomposition MWD 7.24 5.25 5.29 5.93 pyramid. This will not be the case in general; thus, the reaIHS 8.46 7.75 9.22 8.48 son for Step l. Another approach would be to decompose the two images to equivalent spatial resolution levels and then resize the image approximation to fit correctly into the other decomposition pyramid. TweDimensional Wavelet Transform Landsat-SPOT Merging The discrete two dimensional wavelet transform ( ~ D W T ) is For our data set, we used a Landsat TM image of Albuquerjust an extension of the one-dimensional case. As in the oneque, New Mexico, acquired on 15 August 1992. The SPOT dimensional case, the original image is reduced in resolution panchromatic data were acquired over the same area on by a low-pass filter and subsampling, but the images November 1993. The two data sets were registered to within are three-fold. he three detail images are a set of independ0.25 pixels RMS using control points and a first-order polynoent, spatially oriented frequency channels that detail vertical mial fit, the ~~~d~~~ data being reSampled. ~h~ satellite data high frequencies, horizontal high frequencies, and cross-diwere roughly corrected for the atmosphere by using rectional high frequencies (Mallat, 1989). known low reflecting materials or shadows in the scene. A linear stretch was then applied to each channel separately to MWD Image Merger fill the full data range. 1024 by 1024 sections of the data The concept of image merger using MWD arose from the use were used in the merging procedures described below. of the wavelet transform in data compression applications We have shown that the final merged image is depend(Antonini et al., 1990; Froment and Mallat, 1992; Lewis and ent on the wavelet basis selected (Yocky, 1995). For the purKnowles, 1992; Rosiene and Greenshields, 1994). Because pose of this paper, we use the family of compact and well wavelets are used in data compression and reconstruction allocalized wavelets presented by Daubechies (1988). These gorithms, it follows that they may be useful in sensor-comwavelets can be described by the weighting coefficients given pressed information problems. The "sensor-compressed" to h(n). As the number of coefficients used for h(n) increases, information problem and how this problem is addressed by the wavelet becomes smoother. We designate the wavelet by wavelet image merging is presented next. the number of coefficients, i.e., DAUB4 wavelet is the DaubeIn remote sensing, due to optical diffraction or signal-tochies wavelet basis with four weights for h(n). noise design limits, one sensor may provide high spatial resolution at the expense of a wide spectral bandwidth and Standard MWD Image Merger another sensor may have higher spectral fidelity at the exThe "standard" MWD image mergers are generated using the pense of spatial resolution. This is the case with the SPOT steps presented above. TM Bands 3, 2 , and 1, and the regispanchromatic (10.0-m resolution) and the Landsat TM (28.5tered SPOT image were merged using the DAUB4 wavelet bato 120.0-m resolution). The information from both sensors is sis. The panchromatic and multispectral images were compressed information because the real world covers all decomposed to five different resolution levels, ending with resolutions and a wide energy spectrum. The information is the image approximations of 32 by 32 pixels. The multispeccompressed at the sensor, due to the sensor's imaging chartral MWD is performed for each TM spectral band. At the end acteristics. The information gathered by SPOT has retained of the forward transform, each 1132-resolution TM spectral spatial information, but the spectral information has been approximation was inserted into the 1132-resolution SPOT compressed in a lossy manner. The information gathered by panchromatic MWD pyramid and the inverse transform was Landsat TM has retained spectral information, but the spatial performed, giving three separate merged bands: red, green, information has been compressed; again, a lossy compresand blue (RGB). A linear stretch was then applied to the imsion. Yet the information compression for the sensor pair is ages. in different information bands; spatial and spectral. Because The MwD merger is compared with the IHS merger of the lossy compression of these information bands, the in(Haydn et al., 1982). In the IHS merger, the TM bands are dividual decompression of each sensor's information to a transformed from RGB into intensity, hue, and saturation higher resolution level is not possible. However, using the components (Smith, 1978). In the IHS space, the high resoludiffering high resolution information bands (i.e., spatial, tion, panchromatic image is substituted as the intensity comspectral), a mutual (merged) decompression may be possible, ponent. The image is then transformed back into the RGB because information compressed in one sensor is preserved space, thus providing a merger between the multispectral at a higher resolution in the other sensor. and the panchromatic images. The "standard" M ~ D merging (Yocky, 1995) of two data Plate 1 shows the outcome of both the IHS merger and sets can be accomplished by the following steps: the MWD merger. Plate l a is the original TM image, Plate l b (1) The two original images must be spatially registered. Originally, may not be the same array size, so make them TABLE 2. CORRELATION COMPARISON BETWEEN THE ORIGINAL TM, SPOT, MWD the same size (dimensionally square, each side being a MERGED, AND IHS MERGED IMAGES. power of two) by interpolation or replication of pixel values. Correlation (2) Choose the wavelet basis for the transform and choose the final resolution for the MWD. The final resolution should be Image Red Green Blue the same for each MWD pyramid. TMIPAN 0.806 0.805 0.778 (3) Perform the MWD on both images. MWDIPAN 0.951 0.949 0.939 (4) Extract the desired sensor image approximation from its deIHS/PAN 0.981 0.988 0.974 composition pyramid and totally replace the approximation MWD/TM 0.864 0.869 0.848 image in the other sensor's decomposition pyramid. IHSITM 0.809 0.827 0.797 (5) Perform the inverse MWD on the image combination. 1070 September 1996 PE&RS TABLE 3. CORRELATION COMPARISON BETWEEN ORIGINAL TM, SPOT, MWD MERGED, IHS MERGED, AMWD MERGED, AND SELECTIVE RESOLUTION MWD MERGERS.